%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calibrate Risk-Aversion Parameter 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% compute standard errors using bootstrap

%%% 07/20/20 updates: 
% use real accidents
% use the full sample to estimate (Y0, beta)
% (1) keep delta, estimate a model with competeing sources
% (2) set delta=1, estimate change of risk-aversion




clear
clc

path='/Users/.../matlab/code_072020';

result=strcat(path,filesep,'result_072020');
mkdir(result);


load([path,filesep,'bootstrap_data_072020.mat']);



%% Set up 


options=optimset('MaxFunEvals', 1e4, 'MaxIter', 1e4, 'TolX', 1e-8, 'TolFun', 1e-8, 'Display', 'none');   
options_y=optimset('MaxFunEvals', 1e4, 'MaxIter', 1e4, 'TolX', 1e-8, 'TolFun', 1e-8, 'Display', 'none');  

kappa_vec=[1065, 1065*0.2, 1065*0.1];
y_init=zeros(1, 3);

para_save=zeros(size(Y_vec, 1), 7*size(kappa_vec, 2));
obj_save=zeros(size(Y_vec, 1), size(kappa_vec, 2));


para_save2=zeros(size(Y_vec, 1), 6*size(kappa_vec, 2));
obj_save2=zeros(size(Y_vec, 1), size(kappa_vec, 2));


%% Start doing estimation for all subsamples, estimate delta jointly


for rep=1:size(Y_vec, 1)
    
    tic
    
    values_hb1=[Y_vec(rep,:); Y_vec(rep,:)+hb1_coef(rep,:)];
    values_hb2=[Y_vec(rep,:); Y_vec(rep,:)+hb2_coef(rep,:)];

    values_hb1(values_hb1<0)=0.0001;
    values_hb2(values_hb2<0)=0.0001;
    
    
    X=X_vec(:, rep);

    rb_vec_hb1=values_hb1';
    rb_vec_hb2=values_hb2';

    X_new_hb1=[rb_vec_hb1; repmat(X,[1,2])];
    X_new_hb2=[rb_vec_hb2; repmat(X,[1,2])];


    prob_hb1=exp(sum(repmat(logit_coef,[1,2]).*X_new_hb1,1))./(1+exp(sum(repmat(logit_coef,[1,2]).*X_new_hb1,1)));
    prob_hb2=exp(sum(repmat(logit_coef,[1,2]).*X_new_hb2,1))./(1+exp(sum(repmat(logit_coef,[1,2]).*X_new_hb2,1)));


    a_vec=repmat(logit_coef(1:3)',[2,1]); % linear term

    mar_eff_hb1=a_vec.*repmat(prob_hb1',[1,3]).*repmat((1-prob_hb1)',[1,3]);
    mar_eff_hb2=a_vec.*repmat(prob_hb2',[1,3]).*repmat((1-prob_hb2)',[1,3]);


%% Jointly estimate parameters


    for i=1:size(kappa_vec,2)
        
        [rep, i]
        
        kappa=kappa_vec(i);

        values_obs=[values_hb1;values_hb2(2,:)];
        mar_eff=mar_eff_hb1;

        alpha_1=400;

        zeta_init=[1/alpha_1, mar_eff_hb1(1,2)/mar_eff_hb1(1,1)*(values_hb1(1,2)/values_hb1(1,1))/alpha_1, mar_eff_hb1(1,3)/mar_eff_hb1(1,1)*(values_hb1(1,3)/values_hb1(1,1))/alpha_1];

        para_init=[log(0.1*1e-2), -2, -1, -2.5, log(zeta_init)];       


        [para_temp, obj_temp]=fminsearch(@(para) obj_para_all_delta(para, kappa, values_obs, logit_coef, X, y_init, options_y), para_init, options);

        rho=exp(para_temp(1));
        change=exp(para_temp(2))/(1+exp(para_temp(2)));
        change_2=exp(para_temp(3))/(1+exp(para_temp(3)));
        delta=1+exp(para_temp(4))/(1+exp(para_temp(4)))*2;
        zeta=exp(para_temp(5:7)); 

        para_save(rep,(i-1)*7+1:i*7)=[rho, change, change_2, delta, zeta];
        obj_save(rep, i)=obj_temp;
    end
    toc
    
end


para_est=para_save(1,:);
bootstrap_std=std(para_save(2:end,:),0, 1);

save([result,filesep,'results_0720_boostrap1.mat'],'para_save', 'obj_save', 'para_est', 'bootstrap_std');




%% Start doing estimation for all subsamples, fix delta

delta=1;

for rep=1:size(Y_vec, 1)
   
    tic
    values_hb1=[Y_vec(rep,:); Y_vec(rep,:)+hb1_coef(rep,:)];
    values_hb2=[Y_vec(rep,:); Y_vec(rep,:)+hb2_coef(rep,:)];

    values_hb1(values_hb1<0)=0.0001;
    values_hb2(values_hb2<0)=0.0001;
    
    
    X=X_vec(:, rep);

    rb_vec_hb1=values_hb1';
    rb_vec_hb2=values_hb2';

    X_new_hb1=[rb_vec_hb1; repmat(X,[1,2])];
    X_new_hb2=[rb_vec_hb2; repmat(X,[1,2])];


    prob_hb1=exp(sum(repmat(logit_coef,[1,2]).*X_new_hb1,1))./(1+exp(sum(repmat(logit_coef,[1,2]).*X_new_hb1,1)));
    prob_hb2=exp(sum(repmat(logit_coef,[1,2]).*X_new_hb2,1))./(1+exp(sum(repmat(logit_coef,[1,2]).*X_new_hb2,1)));


    a_vec=repmat(logit_coef(1:3)',[2,1]); % linear term

    mar_eff_hb1=a_vec.*repmat(prob_hb1',[1,3]).*repmat((1-prob_hb1)',[1,3]);
    mar_eff_hb2=a_vec.*repmat(prob_hb2',[1,3]).*repmat((1-prob_hb2)',[1,3]);


%% Jointly estimate parameters


    for i=1:size(kappa_vec,2)
        
        [rep, i]
        
        kappa=kappa_vec(i);

        values_obs=[values_hb1;values_hb2(2,:)];
        mar_eff=mar_eff_hb1;

        % try different initial values alpha_1=300(50)700
        % use the one that gives the lowest obj
        
        alpha_1=300;    
        if i==3
            alpha_1=450;
        end
        
        
        zeta_init=[1/alpha_1, mar_eff_hb1(1,2)/mar_eff_hb1(1,1)*(values_hb1(1,2)/values_hb1(1,1))/alpha_1, mar_eff_hb1(1,3)/mar_eff_hb1(1,1)*(values_hb1(1,3)/values_hb1(1,1))/alpha_1];
        para_init=[log(0.1*1e-2), -2, -1, log(zeta_init)];       
        [para_temp, obj_temp]=fminsearch(@(para) obj_para_all_delta2(para, delta, kappa, values_obs, logit_coef, X, y_init, options_y), para_init, options);


        rho=exp(para_temp(1));
        change=exp(para_temp(2))/(1+exp(para_temp(2)));
        change_2=exp(para_temp(3))/(1+exp(para_temp(3)));
        zeta=exp(para_temp(4:6)); 

        para_save2(rep,(i-1)*6+1:i*6)=[rho, change, change_2, zeta];
        obj_save2(rep, i)=obj_temp;
    end

    toc
end


para_est2=para_save2(1,:);
bootstrap_std2=std(para_save2(2:end,:),0, 1);

save([result,filesep,'results_0720_boostrap2.mat'],'para_save2', 'obj_save2', 'para_est2', 'bootstrap_std2');















